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We study how different types of blow-ups can occur in systems of hyperbolic evolution equations 
of the type found in general relativity. In particular, we discuss two independent criteria that 
can be used to determine when such blow-ups can be expected. One criteria is related with the 
so-called geometric blow-up leading to gradient catastrophes, while the other is based upon the 
ODE-mechanism leading to blow-ups within finite time. We show how both mechanisms work 
in the case of a simple one-dimensional wave equation with a dynamic wave speed and sources, 
and later explore how those blow-ups can appear in one-dimensional numerical relativity. In the 
latter case we recover the well known "gauge shocks" associated with Bona-Masso type slicing 
conditions. However, a crucial result of this study has been the identification of a second family of 
blow-ups associated with the way in which the constraints have been used to construct a hyperbolic 
formulation. We call these blow-ups "constraint shocks" and show that they are formulation specific, 
and that choices can be made to eliminate them or at least make them less severe. 
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I. INTRODUCTION 

When studying the Cauchy problem of field theories 
in physics, one has to worry about the existence and 
uniqueness of solutions to the system of evolution equa- 
tions considered. In mathematical terms, one is inter- 
ested in working with problems that are well-posed, by 
which one understands that a unique solution exists (at 
least locally), and solutions are stable in the sense that 
small changes in the initial data produce small changes in 
the solution. In this respect, one usually looks for either 
symmetric or strongly hyperbolic systems of equations 
since Cauchy problems for such systems are known to be 
well posed under very general conditions 

In the case of general relativity, the Cauchy prob- 
lem was studied since the 50 's with the pioneering work 
of Choquet-Bruhat 0, and by the mid 80's a num- 
ber of hyperbolic reductions were known (see for exam- 
ple 0. [j. 0- [(| ■ and more recent reviews in 0,0]). Still, 
those reductions played a minor role in numerical rela- 
tivity, where practically all work using the Cauchy ap- 
proach was based on the Arnowitt-Deser-Misner (ADM) 
system of equations |9|, [lfj . Interest in hyperbolic formu- 
lations in numerical relativity started in the early 90 's, 
with the work of Bona and Masso 0, 0, ^| , but con- 
tinued as a small side branch for a number of years. This 
situation remained until Baumgarte and Shapiro showed 
in 01 that a reformulation of the ADM equations orig- 
inally proposed by Nakamura, Oohara and Kojima jl5T ]. 
and Shibata and Nakamura |l6(, had far superior nu- 
merical stability properties than ADM. Baumgarte and 
Shapiro attributed this to the fact that the new formu- 
lation, which has since become known as BSSN, had a 
"more hyperbolic flavor" . This rather informal state- 
ment was later put on firmer ground in [TtI Il8l Il9j | , and 
today it is understood that ADM is only weakly hyper- 



bolic (and thus not well posed) |2C 



whereas BSSN is 



t well poi 

strongly hyperbolic [lg, H3 • 

The recognition by the numerical relativity community 
of the fact that well-posedness is a crucial ingredient for 
having long term stable and well behaved numerical sim- 
ulations (see [2]]]) has led, in recent years, to an explosion 
in the search for ever more general hyperbolic reductions 
of the Einstein evolution equations. At this time many 
such hyperbolic formulations exist, several of which have 
dozens of free parameters (see for example H3| ) • The 
large number of ways in which one can construct strongly 
or even symmetric hyperbolic formulations has taken us 
to a situation where there are now many more proposed 
formulations than numerical groups capable of testing 
them. At the same time, there is a growing realization 
that in some respects well-posedness is not enough, as 
empirically some hyperbolic formulations have proven to 
be far more robust than others. Some work has been done 
on the analytic side trying to understand what makes 
some hyperbolic formulations better suited for numerical 
work. In particular one can mention the work of Shinkai 
and Yoneda |24|, where the propagation of constraints 
for different formulations is studied by linearizing the 
evolution system around the Schwarzschild background 
and looking for the eigenvalues of the evolution matrix in 
Fourier space, and the work of Lindblom and Scheel [25| . 
where the rate of growth of the constraint violation is an- 
alyzed for a family of symmetric hyperbolic formulations 
using the fact that for such systems one can construct an 
"energy norm". The consensus is that one should look 
beyond the principal part of the system, and study the 
effect of the source terms on the stability. 

In this paper, we want to focus on a different aspect 
that can differentiate between hyperbolic formulations 
and that has been so far overlooked. Well-posed formu- 
lations are known to have well behaved solutions locally, 
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but there is no guarantee that these solutions can exist 
beyond a certain finite time. In fact, on physical grounds 
we expect solutions to fail after a finite time in some 
circumstances, due for example to the formation of sin- 
gularities in gravitational collapse. But there is another 
way in which solutions can become singular after a finite 
time, the best example of which is the formation of shock 
waves in hydrodynamics. In general relativity, and par- 
ticularly in vacuum, we do not expect these type of shock 
waves to develop. Nevertheless, one must remember that 
the evolution equations evolve more than just the physi- 
cal degrees of freedom. In particular, there are also gauge 
degrees of freedom that can cause coordinate singularities 
to arise during the evolution. In [2f| one of us showed 
that coordinate singularities caused by the crossing of 
the characteristic lines associated with the propagation 
of the gauge can in fact easily form. These so-called 
"gauge shocks" are now known to be responsible for the 
fact that some gauge choices are much better behaved 
than others. In particular, in reference |2]J it was shown 
that the well known "1+log" slicing condition, which em- 
pirically has been found to be very robust, is in fact the 
only member of its family that avoids gauge shocks ap- 
proximately. More recently it has been found that, for 
the evolution of Brill wave spacetimes that are very close 
to the critical threshold for black hole formation, the use 
of shock avoiding slicing conditions is crucial 42] . 

But there are degrees of freedom other than the phys- 
ical and gauge modes that also appear in the evolution 
equations of general relativity. These extra degrees of 
freedom have to do with the violation of the constraints, 
and even though for physical initial data they should van- 
ish identically, truncation errors make their presence un- 
avoidable in numerical simulations. It is therefore very 
important to understand how these constraint violating 
modes behave analytically. A main result of this paper is 
the recognition that constraint modes can also give rise to 
the development of blow-ups in a finite time that are very 
similar to the gauge shocks studied before. These blow- 
ups are a property of the specific form of the evolution 
equations and their effects can be significantly reduced if 
one chooses carefully how the constraints are used when 
constructing a hyperbolic system. In this paper we show 
how blow-ups can arise in spherically symmetric relativ- 
ity, and how they can best be avoided by modifying the 
evolution system. We believe that the study of such con- 
straint shocks might help us understand why otherwise 
well-posed and "nice" formulations might behave poorly 
in numerical simulations. 

A comment about our terminology is in order. Bor- 
rowing the language from hydro-dynamics, throughout 
the paper we will refer in a somewhat loose way to blow- 
ups as "shocks" , though this term strictly only refers to 
blow-ups caused by the crossing of characteristic lines. 

This paper is organized as follows. In Section [B] we in- 
troduce the concept of hyperbolicity and describe two dif- 
ferent criteria that can be used to determine when blow- 
ups in the solutions to hyperbolic systems of equations 



can be expected. Section IIIII introduces the simple one- 
dimensional wave equation with sources and a dynamic 
wave speed as an example of how these blow-ups develop. 
In Section IIVI we apply the blow-up criteria to the evo- 
lution equations of general relativity. We start with the 
simple case of 1+1 dimensions, where we recover the well 
known gauge shocks, and later study the case of spherical 
symmetry where we find that constraint shocks can also 
arise. We conclude in Section IVl 



II. HYPERBOLICITY AND SHOCKS 

The system of evolution equations we are interested 
in analyzing are the evolution equations for the Cauchy 
problem of general relativity. In particular we are inter- 
ested in studying the appearance of singular non-physical 
solutions. Such an analysis can be best made using the 
characteristic structure of hyperbolic systems, so we will 
start from the definition of hyperbolicity. We will also 
concentrate on systems with only one spatial dimension 
as this makes the analysis so much simpler. The im- 
portant point of what happens in the multi-dimensional 
case is a matter for future research. Notice that one- 
dimensional systems are in fact relevant in general rel- 
ativity, as they can represent the evolution of systems 
with, for example, spherical symmetry. 

There is one important point that should be men- 
tioned. Throughout this section, and in the rest of the 
paper, we manipulate differential equations by assum- 
ing that partial derivatives in time and space commute. 
This type of manipulation leaves smooth ( "classical" ) so- 
lutions unchanged, but can easily change the speed of 
propagation of shock waves • Still, in this paper we 
will only be interested in smooth solutions, and we will 
consider the development of a shock as a pathology. Our 
whole emphasis is in finding ways to avoid shocks. 



A. Hyperbolic systems 

We will consider quasi-linear systems of evolution 
equations that can be split into two subsystems with the 
following structure 

d t u = -Mid x u-M 2 u K, (2.1) 
d t K + d 2 x u + M 2 K d x K = p K {u, d x u, K) . (2.2) 

Here u and K are n and m dimensional vectors respec- 
tively, and M^' 2 and M^ 2 are matrices whose coefficients 
may depend on the it's, but not on the K J s. 

In order to have a first order system we will intro- 
duce the spatial derivatives D t := d x Ui as extra indepen- 
dent variables, whose evolution equations are obtained 
directly from those of the u's, 

d t D + 8 X (Ml D + M 2 u K)=Q. (2.3) 
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In the following we will always assume that the initial 
data satisfies the constraints Di — d x Ui. This implies 
that spatial derivatives of the u's can always be substi- 
tuted for D's and treated as source terms. 

Let us now define the (n + to) dimensional vector 
v := (D, K). We can then rewrite the system of evolution 
equations as 



d t u = -M(u) v , 

d t v + A(u) d x v = q v {u,v) 



(2.4) 
(2.5) 



where M and A are n x (n + to) and (n + m) x (n + m) 
matrices, 

(Ml M 2 U \ 
M=(M;M5), A^ MiM J, (2,, 

and where the source vector q v is given by 



n 

£A [(du, 



ml)D + {d u Ml)K 



i=l 



,Pk 



(2.7) 



In our primary example, the Einstein equations, the 
vector u consists of gauge variables and components of 
the 3-metric, whereas v contains both variables associ- 
ated with the spatial derivatives of the gauge variables 
and metric components (the D's) and also variables aris- 
ing from the time derivatives of the metric components 
(the K's). Note, furthermore, that the source terms q v 
appearing on the right-hand side of (|2.5|) are in general 
functions of both the it's and u's (typically quadratic on 
the v's). 

The system of equations above will be hyperbolic if the 
matrix A has (n + m) real eigenvalues A^. Furthermore, 
it will be strongly hyperbolic if it has a complete set of 
eigenvectors £j, 

A Ci = AiS • (2.8) 
If we denote the matrix of column eigenvectors by R, 



R — (£l ' ' ' £(n+m) J ) (2-9) 

then the matrix A can be diagonalized as 

R X AR = diag [X lr ■ ■ , A (n+m) ] = A . (2.10) 

Notice that for systems with only one spatial dimension 
the otherwise important distinction between strongly and 
symmetric hyperbolic systems does not arise. 

For a strongly hyperbolic system we define the eigen- 
fields as 

w = BT 1 0. (2.11) 

By analyzing the time evolution of the eigenfields we now 
want to study by which mechanisms blow-ups (i.e. sin- 
gularities) in the solutions can arise. As pointed out 



in [28(, there are basically two different blow-up mech- 
anisms which are, somewhat misleadingly, referred to as 
"geometric blow-up" and the "ODE-mechanism" . Since 
in the first case the derivative of an evolution variable, 
and in the second case an evolution variable itself, be- 
comes infinite within a finite time, the names "gradient 
catastrophe" [2|j and "blow-up within finite time" are 
probably more appropriate. In the following sections we 
will explain the basic idea behind these two mechanisms 
using as a prototype a simple scalar equation, and we will 
show how both mechanisms arc in fact closely related. 
We will also point out how these mechanisms generalize 
to systems of PDE's. 



B. Geometric blow-up and linear degeneracy 

1. Scalar conservation laws 

The first mechanism responsible for blow-ups involves 
only quasi-linear systems of equations. Here the solu- 
tion u under consideration has a well-defined limit at a 
given point and only the derivatives of u become infinite 
there. Typical examples of this situation are obtained 
when solving scalar conservation equations of the form 



d t u + d x F(u) — d t u + a(u) d x u = , 



(2.12) 



where a(u) := dF(u)/du. The blow-up is then due to the 
focusing of characteristics at a point, and the mechanism 
is referred to as "geometric blow-up" . Taking the sim- 
plest nonlinear function F(u) = u 2 /2, we obtain Burgers' 
equation 



dtu + u d x u = , 



(2.13) 



which is frequently discussed in the literature as an exam- 
ple of a genuinely nonlinear PDE leading to shock forma- 
tion (see for example [Hill). The solution of Eq. (|!H3|) 
can be interpreted as a time-dependent one-dimensional 
velocity field u. The equation then states that the char- 
acteristics (i.e. the "flow lines") have zero acceleration, 
that is du/dt — d t u + (dx/dt) d x u = dtu + u d x u = 0, 
which means that particles following those trajectories 
move with constant velocity u = dx/dt. However, unless 
the initial velocity distribution uq(x) is a non-decreasing 
function of x (so that the particles "spread out"), even- 
tually a particle with higher velocity will collide with one 
ahead of it having a lower velocity. In particular, as par- 
ticles initially at rest are not moving at all, u is forced 
to become singular at a finite time if its initial velocity 
distribution has compact support (except in the trivial 
case when uq(x) vanishes everywhere). 



2. Linear degeneracy 

The previous argument coming from Burgers' equa- 
tion H2.13f) can be easily generalized to the case of 
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Eq. H2.12fl . Consider two locations x\ and a; 2 with 
X\ < X2 and corresponding initial values ui — uq{x\) 
and it2 = uo(x2). Just as before, the values of u are con- 
served along characteristics, so unless a(u(x)) is purely 
increasing in x, it is always possible to find locations such 
that a(ui) > 0,(112), so the characteristic lines on the left 
go faster than those on the right (as a(iii) represents their 
velocity). One may readily verify that the lines will in- 
tersect at a time given by 



t* 



X\ - x 2 



a(m) — a(ii 2 ) 



(2.14) 



At the point of intersection u has to take both values u\ 
and tt2, so a unique solution ceases to exist. When this 
happens the spatial derivative of u becomes infinite and 
the differential equation breaks down, in other words no 
smooth solution of (|2.12|) exists after t — t* . 

For smooth initial data, taking the limit \xi — x^ \ — > 
in the expression for t*, we see that this gradient catas- 
trophe will occur at a finite time 



t* = - 



1 



min [d x a(u )] min[a'{u Q )d x u ] 



(2.15) 



in the "future" if initially we have d x a(uo{x)) < 0, 
whereas the problem arises in the "past" if 
d x a(uo(x)) > holds initially. Since in general one 
can not guarantee that every initial data set one would 
like to use satisfies such a condition, the criteria for not 
forming a shock demands that the function a(u) should 
be linear, that is a'(u) = 0. 

The above argument can also be generalized to 
(strongly) hyperbolic systems of equations, see j^J. For 
such systems, the condition one needs in order to avoid 
the formation of shocks associated with the propagation 
of a given eigenfield Wi is for this eigenfield to be "linearly 
degenerate" , which means that the eigenvalue Xi associ- 
ated with Wi must be constant along integral curves of 
the corresponding eigenvector 



dwi 



(n+m) 

EOAi OVj 



dvj dwi 



= 



(2.16) 



C. ODE-mechanism and the source criteria 

1. ODE's with quadratic sources 

In ODE's, PDE's and systems of PDE's, an evolution 
variable itself can become infinite at a point by a process 
of "self-increase" in the domain of influence leading to 
this point. In a somewhat misleading way, the underlying 
mechanism has been given the name "ODE-mechanism" 
(see |28(), since prototype examples are based on simple 
ODE's such as 



dii 
~dt 



cu 



c = constant ^ 



(2.17) 



For non-trivial initial data the solution of <|2.17l) is 



u(t) 



u 



UQ^O 



(2.18) 



1 — UQCt 

This solution clearly blows up at a finite time given by 



t* 



1 

mo c 



(2.19) 



either in the past or in the future, depending on the sign 
of uqc. One can also expect such blow-ups to happen in 
the case when c is not a constant but instead a function 
of time. If c(t) is bounded, one can apply theorems 1 
and 2 of Supposing that the function c(t) satisfies 
the inequality < C < c(t) for < t < T, and that uq 
is positive, then T < 1/(uqC) since u(t) for all positive 
t is bounded from below by uo/(l — uoCt). Similarly, 
supposing that e(t) satisfies the inequality \c(t)\ < C, 
then the initial value problem has a solution for at least 
\t\ < l/\u C\. 



2. The source criteria for avoiding shocks 

Let us now go back to our original system of equations 
(E3|)-(E3. Multiplying Eq. {2} from the left by Rr 1 
we find 

d t (R _1 w) + (R^AR) d x (R _1 w) 

= R-y*- + [dtR.- 1 + (R _1 AR) ^R" 1 ] v , (2.20) 

which, by making use of H2.10|) and (|2.11|) . yields 

d t w + Ad x w = q w , (2.21) 

where 



R 



+ [3 t R _1 + Ac^R -1 ] 



(2.22) 



In this way we have obtained an evolution system 
where on the left-hand side of H2.21J1 the different eigen- 
fields Wi are decoupled. However, in general the equa- 
tions are still coupled through the source terms q Wi . In 
particular, if the original sources were quadratic in the 
i>'s, we will have 



dw 

~dt 



— = d t Wi + Kd x Wi 



(n+m) 

cijkWjWk - 

3,k=l 



■0{w) , (2.23) 



where d/dt := dt + Xid x is the derivative along the corre- 
sponding characteristic. As pointed out for a similar sys- 
tem in |38l l34j| . here the Cmw\ component of the source 
term can be expected to dominate, so mixed and lower 
order terms can be neglected. Though we have no proof 
of this statement in the general case, one can expect it to 
be true at least for systems with distinct eigenspeeds, as 
mixed terms will then be suppressed when pulses moving 
at different speeds separate from each other. The effect 
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of the term cmwf, on the other hand, will remain even 
as the pulse moves. The numerical simulations presented 
in the following sections show empirical evidence that 
reinforces this argument. 

We can then rewrite the above equation as 



dwj 
~dt 



(2.24) 



which has precisely the form of the ODE studied in the 
previous section. In order to avoid a blow-up one would 
then have to demand that the coefficients cm vanish. We 
call this the "source criteria" for avoiding blow-ups. 

There is a very important property of our system of 
equations regarding the coefficients cm that come from 
the source terms q w . From H2.22(l one could expect con- 
tributions to these coefficients coming both from the orig- 
inal sources q v and from the term in brackets involving 
derivatives of R _1 . However, one can show that this is 
not the case and for the systems under study the contri- 
butions to cm coming from the term in brackets cancel 
out, that is, all contributions to cm come only from the 
original sources q v . 

In order to see this we start by rewriting the term 
inside brackets on the right hand side of l|2.22|l as 

n 

dtR-i+AdeBT 1 =J2{d t u l +Ad x u l )d Ul IL- 1 . (2.25) 
i=i 

Since both the time and space derivatives of ui can be 
written in terms of w's, and the above term multiplies 
the vector v in Eq. I|2.22|l . it is clear that this term will 
give rise to quadratic terms in the w's, and hence in the 
w's. The question is whether these quadratic terms will 
produce a contribution to the coefficients cm. From the 
last equation it is clear that no such contribution will 
exist if the following condition is satisfied 



d 
dwj 



(dm + Xid x ui) = 0, V i < (n + m), I < n. (2.26) 



We will now show that this condition is indeed always 
satisfied. Notice first that, from the definition of the 
eigenfields w and the matrix R, one can easily see that 



d 

dwi 



= & ■ V, , 



(2.27) 



with the eigenvector corresponding to the eigenvalue 
A,;. Now, from equation (|2.4Jl and the definition of the 
D's we have 

(n+m) 

dm + A, d x ui = - M i3 v i + x i D i . ( 2 - 28 ) 



which implies that 



(n+m) 



£• • v, (dm + Xi d x u{) = Xiin - ( 2 - 29 ) 



where £jj is the j component of the vector and where 
we have used the fact that the first n components of v are 
precisely the D's (remember that by construction I < n). 

To finish the proof we now use the fact that the first n 
rows of the matrix A are given by the matrix M, and also 
the fact that £j is an eigenvector of A with eigenvalue Xi , 
which implies 

(n+m) (n+m) 

E M ii&j= E -M"-.. '-V6/. (2.30) 

3=1 3 = 1 



from which we finally find 



& • V„ (dm + A 4 d x u{) = 



(2.31) 



This completes the proof that condition H2.26(l always 
holds (as long as the constraints Di = d x Ui are satisfied), 
which in turn means that the term in square brackets 
in 12.22fl does not contribute to the coefficients cm- 

We want to make another important comment here: 
As the eigenvectors diagonalizing the matrix A are ob- 
tained only up to an arbitrary rescaling, also the eigen- 
fields Wi are not unique. In particular, any Wi can be 
multiplied by an arbitrary function of the it's to obtain 
ibi — fli(u) Wi. However, since the w's are related to 
derivatives of the u's, such a rescaling will introduce new 
quadratic source terms, so one would in general not ex- 
pect the coefficients cm to be invariant under rescalings 
of the eigenfields. 

Remarkably, for the systems of the type (|2.4|I - H2.5|I that 
we are interested in, it turns out that such rescalings 
of the eigenfields have no effect on the coefficients cm. 
The proof of this is again related to condition l|2.26|l . In 
general, if we rescale the eigenfunctions as Wi = ili(u) Wi, 
we will find that 

d t Wi + X i d x w l = Q (d t Wi + Xid x Wi) 

n 

+Wi 2J d Ul fl (dm + X l d x ui) 
i=i 

n 

= Q q Wi + Wi J2d ui n (d t ui + Xid x ui) . (2.32) 
i=i 

From this we see that, although the rescaling does intro- 
duce new quadratic terms, condition l|2.26|) guarantees 
that no new contributions to the coefficient cm will arise, 
i.e. the source criteria for avoiding blow-ups is invariant 
with respect to rescalings of the eigenfunctions (again, as 
long as the constraints Di = d x Ui are satisfied). 



3. Is the source criteria necessary and sufficient in order to 
avoid blow-ups? 

A question one might immediately ask is whether the 
source criteria introduced above is necessary and suffi- 
cient to avoid blow-ups. Although we have no proof at 
this time, numerical experiments (such as those shown in 
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later sections) indicate that whenever the source criteria 
is not satisfied blow-ups do develop, which would sup- 
port our conjecture that the criteria is indeed necessary 
in order to avoid blow-ups. 

About it being sufficient, it clearly is not. This can be 
seen from the following example. Consider the following 
system of two equations, 



d t vi 
d t v 2 



d x v 2 
d x vi 



-2v\V 2 , 



which can easily be diagonalized to find 

d t w+ + d x w + = w 2 , , 
d t w- - d x w- = w% , 



(2.33) 
(2.34) 

(2.35) 
(2.36) 



with w± := v± ± v 2 . The diagonalized system clearly 
satisfies the source criteria. Now, take initial data such 
that v\ = k = const and v 2 — in a large spatial region. 
This implies that in that region w\ — w 2 — k. Since the 
spatial derivatives vanish and both fields are in fact equal, 
the equations above are of the form (|2.17l) and the fields 
will blow up in finite time (provided the region where the 
fields where initially equal is large enough). This initial 
data is clearly very special, but it does show that the 
source criteria is not sufficient in order to avoid a blow- 
up. Nevertheless, for more generic data this situation 
will be very rare. 

We have in fact performed numerical experiments with 
systems of the above form, but generalizing the source 
terms to 



d t w± ± d x w± 



a±w± + b±w + w- + c±wZ. . (2.37) 



When using Gaussians with a small amplitude of order 
0(e) as initial data for w±, then whenever the coeffi- 
cients a± are non-zero one typically finds that blow-ups 
occur on a timescale of order 0(l/e). If, on the other 
hand, a± = and one has only mixed terms and/or terms 
quadratic in the other eigenfield in the sources, blow-ups 
again eventually develop, but now on a timescale of order 
0(l/e 2 ). Hence, if one is interested in propagating small 
perturbations, then satisfying the source criteria should 
allow one to obtain longer evolutions. 



D. Relationship between the different blow-up 
mechanisms 



In order to understand the relationship between the ge- 
ometric blow-up and the ODE-mechanism, we will study 
a system of two variables constructed from the simple 
scalar conservation law l|2.12l) by introducing either the 
time or space derivative of the function u as an extra 
independent variable. 

We start by introducing D := d x u as a new variable. 
One then obtains the system 

d t u = -a(u) D , (2.38) 
d t D + a(u) d x D = -a!{u) D 2 , (2.39) 



where the evolution equation for D has been found by 
differentiating l|2.12f) with respect to x and exchanging 
the order of d t and d x . 

As we are interested in studying solutions of the origi- 
nal scalar conservation law, but seen from a different per- 
spective, will only consider initial data such that the con- 
straint D := d x u is satisfied. Remembering that along a 
characteristic line of (|2.12|) u is constant (and hence also 
a(u) and a'(u)), the equation 



dD 
~dt 



d t D + a(u) d x D = -a'(u) D 2 , 



(2.40) 



arising from (|2.39|l can be easily integrated. We find that, 
along the characteristic, the following relation holds 



D(t) 



Do 



1 + D a'(u)t 

This clearly becomes infinite at a time t* given by 



(2.41) 



D a'(u) 



(2.42) 



Let us now introduce K := dtu instead of D as an 
extra variable. We then obtain the system 



d t u = K , 
d t K + a(u) d x K 



a'(u) 
— — K 

a(u) 



(2.43) 
(2.44) 



Here the evolution equation for K has been derived by 
taking a partial derivative with respect to t of (|2.12() . As 
before, by integrating the equation 



dK 



u d t K + a(u)d x K 
at a(u) 

along the characteristic one finds 



a '( u ) R 2 



K(t) 



K 



1 - (K a'(u) /a{u)) t 
which diverges at a time given by 



t* = 



a(u) 
Koa'(u) 



(2.45) 



(2.46) 



(2.47) 



These two examples are nothing more than our origi- 
nal scalar equation 12.12fl in disguise. However, they are 
in fact linearly degenerate by the definition given above 
as the only eigenvalue a(u) is independent of D and K, 
respectively. They still give rise to a blow-up, as they 
should, but this time the blow-up appears through the 
ODE-mechanism instead of the original geometric blow- 
up mechanism. Notice that, from l|2.42|l and l|2.47|l . one 
can conclude that a condition for not having a blow-up 
in finite time is a'(u) = 0, which is the same condition 
we found in Sec. Ill B 21 above. This shows clearly that 
what can be considered a geometric blow-up of a given 
variable u can always be reinterpreted as an ODE-type 
blow-up of its derivatives, so both blow-up mechanisms 
are closely related. 
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E. Indirect Linear Degeneracy 

Linear degeneracy turns out to be insufficient for avoid- 
ing blow-ups in the particular case of system I|2.4|) - (I2.5I) 
for two reasons. The first reason has to do with the 
presence of non- vanishing source terms q v and has been 
discussed in Sec. Ill (Jl above. The other reason is simply 
the fact that for these type of systems the eigenvalues of 
the characteristic matrix A depend only the u's and not 
on the v's, which means that all eigenfields are linearly 
degenerate in a trivial way. We have already seen an 
example of this in the previous section when we consid- 
ered the simple scalar conservation law and introduced 
derivatives as extra independent variables. 

For this reason the concept of "indirect linear degen- 
eracy" was introduced in |2g. This simply replaces the 
eigenvalue A^ in equation (|2.1(il) by its time derivative: 

• (n+m) • 

. • r* UAj x — > CM, Ova „ 

^•6 = -^ = V j^jt- = 0. 2.48 

This new condition yields non-trivial results for the sys- 
tem l|2.4ll - H2.5|) if the time derivatives of the u's, appear- 
ing when differentiating A^ with respect to time, depend 
on the corresponding 

It is in fact not difficult to see where the indirect linear 
degeneracy condition comes from. Consider the system 
of two equations 

d t u = p(u,v,d x u) , (2.49) 
d t v + X(u) d x v = q(u, v, d x u) , (2.50) 

with p linear in v and d x u. We now extend the above sys- 
tem by introducing the variable D := d x u. This means 
that the sources p and q are now functions of (u,v,D). 
The full system will then be 

d t u = p , (2.51) 
d t D - d x p = 0, (2.52) 
d t v + A(it) d x v — q , (2.53) 

which is exactly of the form 12.4|) - (|2.5|l . Let us for a 
moment assume that q — 0. In that case it is clear 
that v will be constant along the characteristics lines 
x = xq + X(u) t. The simplest example is obtained 
when X(u) — u and p = v — uD, since then we find that 
along the characteristics du/dt = v (provided that the 
constraint D = d x u remains satisfied). This means that 
along those lines we have u — Uq + vt, so the character- 
istics have constant acceleration given by v (since u is 
the characteristic speed). If initially vq(x) = v(t = 0,x) 
has negative slope in a given region, the characteristics 
are then guaranteed to cross (as lines behind accelerate 
faster than those in front). At the point where this hap- 
pens the gradient of v will become infinite and we will 
have a blow-up. For cases when p is a different function 
one can not integrate the equations exactly, but the same 
general idea will hold. Of course, when the source term 



q is not zero one could imagine that q can be chosen in 
such a way as to avoid the crossing of characteristics, 
but such a choice would clearly not be generic. The only 
way to be sure that there will be no blow-up is to ask 
for dp/dv — 0. Indirect linear degeneracy is simply the 
generalization of this condition to the case of a system 
with more equations. 

The argument given above, however, is clearly not rig- 
orous. Indirect linear degeneracy is therefore still a more 
or less ad hoc condition. Part of the reason for discussing 
it here is precisely to study its relevance in different cases 
by numerical experiments. As our results in the following 
sections show, indirect linear degeneracy and the source 
criteria often yield the same conditions for avoiding blow- 
ups. When they do not, the source criteria seems to 
be more important. Exploring the link between indirect 
linear degeneracy and the source criteria is something 
that should be further investigated, and we are currently 
working in that direction. 

III. THE WAVE EQUATION WITH SOURCES 
AND A DYNAMIC WAVE SPEED 

A. Blow-up formation 

As an example for the type of evolution systems stud- 
ied in the previous sections we will consider the simple 
scalar wave equation with sources, 

dfu — c 2 (u)d x u = q(u,d t u,d x u). (3.1) 

Here we allow the wave speed c to be a function of the 
wave function u. The source term q, on the other hand, 
can depend both on u and its first derivatives. Intro- 
ducing D = d x u and K = d t u, we can rewrite the wave 
equation as 

d t u = K , (3.2) 
dtD - d x K = 0, (3.3) 
8 t K - c 2 d x D = q, (3.4) 

which is of the form H2.4fl - H2.5|) . One can readily verify 
that the eigenvalues of the characteristic matrix, 

A^(_° c2 - 1 ), (3.5) 

are A± = ±c, with corresponding eigenfields (the nor- 
malization is chosen for convenience) 

w± = 1 (K T cD) . (3.6) 

The linear degeneracy criteria then is trivially satisfied 
since the eigenvalues depend only on u. However, when 
we calculate the time derivative of the eigenvalues we find 

A± = ±c'd t u = ±c'K = ±d {w+ + w-) . (3.7) 
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The indirect linear degeneracy condition l|2.48[l asks for 
the derivatives d\±/dw± to vanish, which implies 







(3.8) 



This is no longer trivial and corresponds to what one 
would in fact expect: The wave speed should be inde- 
pendent of the wave function if we want no shocks to 
develop. 

Let us now turn to the source criteria. We find 



dw± 1 c' 



W+W- — w n 



(3.9) 



Notice that the evolution equation for a given eigenfield 
w± contains no quadratic terms on itself other than those 
that might come from q (it has only mixed terms and 
terms quadratic in the other eigenfield). The source cri- 
teria then demands that the source term q should be free 
of the quadratic terms and w_. If we assume that q 
is of the form 



q = A c 2 D 2 + B c DK + CK 2 

= {A-B + C)w 2 + -2{A-C)w + w. 
+ {A + B + C)w 2 _, 



(3.10) 



with A, B and C arbitrary functions of u, it follows that 
in order to avoid blow-ups these functions have to satisfy 



B = 



A + C = 



(3.11) 



B. Numerical results 



We have performed a series of numerical simulations 
for wave equations which satisfy or violate indirect lin- 
ear degeneracy and/or the source criteria. The results 
from these simulations are summarized in Table [I] All 
simulations have been performed using a method of lines 
with fourth order Runge-Kutta integration in time, and 
standard second order centered differences in space (with 
no artificial dissipation added |43j). 

As initial data we have taken u(t = 0) = 1 and hence 
D(t = 0) = 0, together with the derivative of a Gaussian 
for the time derivative of u, i.e. 



K(t = 0) = - (2kx/<t 2 ) exp (-x 2 /a 2 ) 



(3.12) 



Here we used the derivative of a Gaussian and not a sim- 
ple Gaussian in order to excite perturbations where u is 
both smaller and larger than its initial value. 

For all the runs shown here we have used the particu- 
lar values k = 0.1 and a = 0.3. These rather strong and 
localized perturbations are motivated by the fact that we 
wanted to see shock formation early, in particular before 
the variable u changes sign (since the first two examples 
with eigenspeeds given by ±it are not strongly hyper- 
bolic for u = 0). However, also for other values of k and 



i.l.d. 


s.c. 


(c, q) 


result 


no 


no 


c — u, q — 2uD 2 


blow-up (Fig.0 


no 


yes 


c = u,q = 2{uD 2 -K 2 /u) 


blow-up? (Fig.|2Jl 


yes 


no 


c = l,q = 2D 2 


blow-up (Fig.[3Jl 


yes 


yes 


c = l,q= 2(D 2 ~ K 2 ) 


no blow-up (Fig. 



TABLE I: Summary of results from simulations of a wave 
equation that satisfies (yes) or violates (no) both indirect lin- 
ear degeneracy (i.l.d.) and the shock criteria (s.c). 



cr we have seen qualitatively very similar behavior (and 
whether or not u crosses zero does not seem to play a 
role in the formation of shocks). 

For the runs with highest resolution we have used 
80, 000 grid points and a resolution of Aa; = 5 x 10~ 4 , 
which places the boundaries at ±20, together with time 
steps of At = Ax/2. In addition, for each evolution vari- 
able we have computed the convergence factor r\ which, 
using three runs with high (u h ), medium (u m ) and low 
(u l ) resolutions differing in each case by a factor of two, 
can be calculated as e.g. 



N 



1_ Y^J 
V, Z^7 = l 



(3.13) 



In the plots we show four different convergence fac- 
tors: We denote with a triangle the convergence factors 
obtained when comparing runs with 80,000, 40,000 and 
20, 000 grid points and a spatial resolution of 5 x 10~ 4 , 
10~ 3 and 2 x 10~ 3 . We use boxes, diamonds and stars to 
denote the convergence factors when gradually lowering 
all three resolutions by a factor of two. For second order 
convergence we expect r\ ~ 4. 

The first test we have done corresponds to a case that 
violates both blow-up criteria. We obtain such a system 
by simply taking a time derivative of Burgers' equation. 

u 2 d 2 u = 2u(d x u) 2 and hence identify 



We find d 2 u 



2uD 2 



(3.14) 



Results for this simulation can be found in Figure ^ In 
the left panel we show snapshots of the evolution of the 
variables u, D and K in steps of At — 1, and in the 
right panel we show convergence factors for the previ- 
ously mentioned series of different resolutions. From the 
figure we clearly see that, as expected, shocks do form, 
with large gradients developing on u and large peaks on 
D and K. Moreover, from the convergence plots we see 
that there is a clear loss of convergence, and as the reso- 
lution is increased, this loss of convergence becomes more 
sharply centered around a specific time t ~ 7, indicating 
that the blow-up happens at this time. 

For the second example we have chosen a situation 
where indirect linear degeneracy is violated but the 
source criteria holds. Numerical results for the case 



q = 2{uD 2 - K 2 /u) 



(3.15) 
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FIG. 1: Results for the case when both criteria are violated. 
This simulation corresponds to the time derivative of Burgers' 
equation, for which we have c = u and q = 2uD 2 . In the left 
panel a sharp gradient is clearly visible on u, and large peaks 
can also be seen on D and K. In the right panel we see a 
clear loss of convergence that becomes sharper as resolution 
is increased in the order "star, diamond, box and triangle", 
indicating a blow-up at a time t ~ 7. 

are shown in Figure[3] (a simpler case arising for a vanish- 
ing source term, q — 0, yields very similar results). The 
figure shows large peaks developing in both D and K, 
and a sharp gradient developing in u. The convergence 
plots show some loss of convergence for the lower reso- 
lutions but, in contrast with the previous example, con- 
vergence seems to improve with resolution. This would 
seem to indicate that although sharp gradients do de- 
velop, a real blow-up has not occurred. Nevertheless, 
such sharp gradients are difficult to resolve numerically, 
so their presence is undesirable. 

The third simulation corresponds to a case that sat- 
isfies indirect linear degeneracy, but violates the source 
criteria, with wave speed and source term given by 

c = l, q = 2D 2 . (3.16) 

Results for this run are shown in Figure |31 As before, 
we see that both D and K are developing large peaks. 
The evolution variable u is developing both a large peak 
and a large gradient. The convergence plots show a loss 
of convergence at lower resolutions that improves as the 
resolution is increased. However in this case all runs crash 
at t ~ 7, indicating that a blow-up has indeed occurred 




FIG. 2: Simulation for a case that violates indirect lin- 
ear degeneracy, but satisfies the source criteria (c = u, 
q = 2{uD 2 - K 2 /u)). The left panel shows that D and K 
develop large peaks, while u develops a sharp gradient but no 
peak. The convergence plots on the right panel indicate also 
loss of convergence in D and K, however convergence seems 
to improve with resolution. 

at that time. 

Finally, our last test corresponds to the case when both 
criteria are satisfied, with the wave speed c and source 
term q given by 

c = l, q = 2(D 2 -K 2 ). (3.17) 

Results for this simulation are shown in Figure^] We see 
that the solution behaves in a wave-like manner, with no 
evidence of a blow-up. This result is reinforced by the 
convergence plots indicating that we have close to second 
order convergence during the whole run for all resolutions 
considered, with no evidence of loss of convergence at any 
time (notice the change in scale with respect to previous 
plots). 

From the previous simulations it is clear that, for the 
scalar wave equation with a dynamic wave speed and 
sources, sharp gradients and blow-ups are only avoided 
when both indirect linear degeneracy and the source cri- 
teria are satisfied. In particular, the case with c = I and 
q = 2(D 2 — K 2 ) behaves very similar to what one would 
expect from the standard wave equation with unit wave 
speed and a vanishing source term. 

This observation can be easily understood by gener- 
alizing an example suggested by L. Nirenberg in |35| : 
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FIG. 3: Results of a simulation for the case that satisfies indi- 
rect linear degeneracy, but violates the source criteria (c = 1, 
q = 2D 2 ). The left panel shows that D and K develop large 
peaks, while u develops both a peak and a large gradient. The 
right panel shows that convergence is lost at low resolutions 
but improves at higher resolutions, until a final time of t ~ 7 
when runs at all resolutions crash, indicating a true blow-up. 

Smooth solutions that extend globally in time will al- 
ways exist if one can find a smooth transformation of the 
form u — u(u), for which u satisfies the standard wave 
equation, dfu — d%u — 0, see [3(||37J. One may readily 
verify that for our particular example with wave speed 
and source term given by l|3.17[) . this is indeed the case 
for the variable u = exp(2u). 



IV. THE EINSTEIN EQUATIONS 

In the previous sections we have described how blow- 
ups can be produced in systems of hyperbolic equations, 
and what conditions need to be satisfied in order for these 
to be avoided. We have also considered one simple exam- 
ple, the wave equation with sources and a dynamic wave 
speed. We will now turn our attention to the system 
we are most interested in, namely the evolution equa- 
tions of general relativity. In this paper, we will restrict 
ourselves to two cases, "toy" 1+1 relativity, and spheri- 
cally symmetric relativity, and leave the important three- 
dimensional case for a future work. 

We believe that it is important to mention here the 
main results which will be presented in this section. 




FIG. 4: Simulation for the case that satisfies both blow- 
up avoiding criteria. In this case we have taken c = 1 and 
q = 2(D 2 -K 2 ), and we find wave-like behavior with no evi- 
dence of a blow-up. 

In the first place, we will recover the results reg arding 
"gauge shocks" discussed already in [2(| |2?], |38| But 
the most important result that we will show is the fact 
that, for the spherically symmetric case, one can also 
identify a second family of blow-ups that are not associ- 
ated with the gauge but rather with the violation of the 
constraints. We will refer to such blow-ups as "constraint 
shocks" , since they are clearly associated with the way in 
which the constraints have been added to the evolution 
equations. These constraint shocks will correspond to 
blow-ups in the hamiltonian and momentum constraints 
at a finite time as the numerical example at the end of 
Sec. II V (J 21 shows gj. 

Since the blow-up analysis assumes that we have a 
strongly hyperbolic system, we will in each case begin 
by constructing such a hyperbolic system for the Ein- 
stein equations. Notice that there is no unique way to 
obtain hyperbolic evolution systems from the Einstein 
equations and we will use this fact to explicitly construct 
formulations that avoid constraint shocks. 



A. Einstein equations in 1+1 dimensions 

Let us assume that we have standard general relativity 
in one spatial dimension. It is well known that in such a 
case the gravitational field is trivial and there are no true 
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dynamics. However, one can still have nontrivial gauge 
dynamics that can be used as a simple example of the 
type of behavior one can expect in the higher dimensional 
case. We will start from the "standard" Arnowitt-Deser- 
Misner (ADM) equations for one spatial dimension [To| . 
where by standard we mean the version of York @ . Since 
we want a hyperbolic system of evolution equations that 
includes the gauge, we will use the Bona-Masso family of 
slicing conditions |Xl| : 



d t a = —a 2 f(a)trK , 



(4.1) 



where / = f(a) > identifies the member of the Bona- 
Masso family being used (for example, / = 1 corresponds 
to harmonic slicing, and / = 2/a to the so-called 1+log 
slicing). For simplicity we will also restrict ourselves to 
the case of vanishing shift vector. 

Following |26| , the two-dimensional vector u will con- 
sist of the lapse function a and the spatial metric function 
.9 := g xx as components. The vector v, on the other hand, 
is a three-dimensional vector with components given by 
the logarithmic spatial derivatives of a and g, and the 
unique component of the extrinsic curvature (with mixed 
indices). That is, 



(D a ,D g , K) 



where 
D n 



dzlna, D g :=d x lng, K := K% 



(4.2) 



(4.3) 



(Note that in |2(j the variable D = d x g/2 is used in- 
stead of D g and K xx is used instead of K%.) The fact 
that we define the D's as logarithmic derivatives instead 
of simple derivatives is in order to simplify the result- 
ing equations, and makes no significant difference in the 
analysis of Sec. [HI 

The ADM evolution equations for the vectors u and v 
turn out to be 



and 



d t a = 
dtg = 

d t D a + 8 X (afK) 
d t D g + d x (2aK) 
d t K + d x (aD a /g) 



-a 2 fK , 
-2agK , 



(4.4) 
(4.5) 



, (4.6) 
, (4.7) 
a (K 2 - D a D g /2g) . (4.8) 

This system has again the form Q2.4[l - I|2.5[l . In par- 
ticular, the last three equations can be written as 
dtv + A(u)d x v = q v , where 



and 



a/ 
2a 
a/g 

-a(f + af)D a K 
-2aD a K 
(K 2 -D 2 Jg + D a Dg/2g) 



(4.9) 



(4.10) 



When studying the characteristic structure of the sys- 
tem of equations above we find the following eigenvalues 



A = , \{ = 

with corresponding eigenfunctions 

/ 




Wo 



w 



f 



D„ 



fgK±D a 



(4.11) 

(4.12) 
(4.13) 



The system is therefore strongly hyperbolic as long as 
/ > 0, with one eigenfield propagating along the time 
lines and the other two propagating with the "gauge 
speeds" A£ = ±ay/ f/g. 

In order to study the possible formation of shocks 
for the propagating eigenfields one can immediately see 
that the direct linear degeneracy criteria as formulated 
in (|2.16|) can not be used, since A£ does not depend on 
either D a , D g or K. The indirect linear degeneracy con- 
dition, however, yields 



dk f ± 



(n+m) 

E 



±- 



a 



1-/ 



of 
2 



0, 



dvj dw f ± 2g 

(4.14) 

where the last step comes from expressing the time 
derivatives of a and g contained in in terms of K 
using (O) and JO). 

For the source criteria, on the other hand, we need to 
determine the term quadratic in w± in the source terms 
associated to the evolution equation for w£ itself. We 
find 



dt 



a 



2Vfg 
o 



i-/ 



of 
2 



/ 2 



f 

WqW ± 



f f 

w ± w T 



(4.15) 



Asking for the coefficient of the quadratic term to be zero 
one finds 



Jff 
c ±±± 



2Vfg 



i-/ 



of 
2 



. 



(4.16) 



It is interesting to note that here both indirect linear 
degeneracy and the source criteria yield the same condi- 
tion for avoiding blow-ups, namely 



J 2 







(4.17) 



The reason why this is so is not completely clear, but it 
probably implies that in this case the sources and char- 
acteristic speeds are not independent of each other. 

The shock avoiding condition l|4.17|l has been studied 
before in [26l |27[ l38| . Its general solution is 



/(«) 



1 



const 



(4.18) 
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Reference ,27] also considers some approximate solutions 
that are more useful for numerical simulations. Notice 
also that, since in this simple case we only have gauge 
dynamics, these shocks are directly associated with the 
foliation itself, and for this reason they are known as 
"gauge shocks" . 



B. Einstein equations in spherical symmetry 

1. Standard ADM equations 

In order to generalize the previous system to spherical 
symmetry, we start with the spatial line element 



dl 2 = A(t, r) dr 2 + r 2 B(t, r) dVl 2 



(4.19) 



where dfl 2 = dd 2 + sin 2 9d4> 2 denotes the usual solid angle 
element. Note that in this case the vector u will consist 
of the lapse a and the metric components A and B. Fur- 
thermore, the w's are given by the spatial derivatives of 
these quantities 

D a = d r In a , D A = d r In A , D B = d r InB , (4.20) 

together with the extrinsic curvature variables (again 
with mixed indices) 



k a = k: 



K B = K% = K* 



(4.21) 



That is, 

u=(a,A,B) , v = (D a ,D A ,D B ,K A ,K B ) . (4.22) 

In the following we will assume for simplicity that we 
are in vacuum and that the shift vector vanishes. Us- 
ing the Bona-Masso slicing condition (|4.1|) , the evolution 
equations for the u's become 



d t a = ~a 2 .f (K A + 2K l 
d t A = -2aAK A , 
d t B = -2aBK B . 



(4.23) 
(4.24) 
(4.25) 



The evolution equation for the u's can be obtained di- 
rectly from the ADM equations and the definition of the 
D's. These equations can again be written in the form 
dtv + A(u)d r v = q v , where the characteristic matrix is 



/ af 2a/ \ 

2a 

2a 

a/A a/A 

V a/2A J 



(4.26) 



and the source terms are given by 

q Da = -a{f + af')D a {K A + 2K B ) 
q Dji = -2aD a K A , 
qn B = -2aD a K B , 



(4.27) 
(4.28) 
(4.29) 



Ik A = —r 



D..(D a -^)-^.(D A -D B ) 



1 



AK A (K A + 2K B )-- (D A - 2D B ) 
r 



(4.30) 



1Kb = -7TT 



a 

■2A 



Dl: ( D a - + D t 



2AK B {K A + 2K b ) + - (2D a - D A + AD B ) 
r 

(4.31) 



2 

r~B 



Furthermore, the Hamiltonian and momentum con- 
straints take the form 



C h 



-drD B + —[D A - — 



AK B {2K A + K B ) 
D A — 3D B A-B 
+ r 2 B 



r 

-d r K B 



0, 



(4.32) 



£f + l)(K A -K B ) = 0. (4.33) 
2 r 



2. ADM equations in new variables 

Rather than working with the standard ADM equa- 
tions described in the last section, following [3^ we will 
introduce the "anti-trace" of the spatial derivatives of 
the metric components, D — D A ~ 2D B , and the trace of 
the extrinsic curvature, K — K A + 2K Bl as fundamental 
variables instead of D A and K A . This choice of variables 
makes the hyperbolicity analysis more transparent. The 
vector v will then be 

v = (D a ,D,D B ,K,K B ) , (4.34) 

and the evolution of the u's will be given by 

d t a = ~a 2 fK , (4.35) 
d t A = ~2aA{K -2K B ) , (4.36) 
d t B = -2aBK B . (4.37) 

For the u's, one again obtains a system of the form 
dtv + A(u)d r v = q v , but this time with characteristic 
matrix 



/ af \ 

2a -8a 

2a , (4.38) 

a/A 2a/ A 

V a/2A 0/ 
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and source terms 

q Da = -a(f + af')D a K 
q D = -2aD a {K - 4K B ) 
Qd b = -2aD a K B , 
a 



qK = ~A 



LI, [D a -j)-D B (D+^ 



-(D a -D + D B )~ — (A-B 



r 2 B 



a 

2A 



D B D, 



D 
2~ 



2AKK B 



- (2D a -D + 2D B ) -^-(A-B) 
r r^B 



(4.39) 
(4.40) 
(4.41) 

- AK 2 
(4.42) 

(4.43) 



Finally, the Hamiltonian and momentum constraints take 
the form 



C h 



d r D 



D B 



2 V 2 J 



AK B (2K - 3K B ) 
D — D B A-B 



Cm 



r 

-d r K B 



r 2 B 

Db 1 
2 r 



, (4.44) 
{K - 3K B ) = . (4.45) 



3. Modifying the equations by using the constraints 

It turns out that, as they stand, neither the matrix A 
of the original ADM system, (|4.26l) . nor the one of the 
rewritten system, ()4.38|l . has a complete set of eigenvec- 
tors for all / > 0, so the systems of evolution equations 
are not strongly hyperbolic. Interestingly, strong hyper- 
bolicity only fails for / = 1, which corresponds to har- 
monic slicing (this only happens in spherical symmetry, 
in the full 3-dimensional case strong hyperbolicity fails 
for ADM much more severely). Since harmonic slicing 
is such an important condition for both theoretical and 
practical reasons, the systems described above are not 
very useful. 

Let us concentrate on the second system of evolution 
equations. By making use of the constraint equations, 
(|4.44|) and (|4.45() . its principal part can be modified to 
construct a strongly hyperbolic system for all / > 0. In 
particular, adding multiples of the constraints will mod- 
ify the third and fifth columns of the matrix l|4.38|) . We 
will consider adjustments to the evolution equations for 
the i;'s of the form 



d t v t 



(n+m) 

Ai-jdrV-j 



+ a (KA n *C h + mi A M *C m ) = q % . (4.46) 



Here the terms {hi,rrii} are allowed to depend on /(a), 
such that for harmonic slicing these coefficients reduce 



to constants. Furthermore, the exponents {Hi,A4i} are 
fixed by looking at the characteristic matrix and demand- 
ing that for / = 1 its entries have homogeneous powers 
of A. Doing this we find 

(4.47) 
(4.48) 
(4.49) 
(4.50) 

and the characteristic matrix then takes the general form 





= 7~Ld = 


Hd b — - 


-1/2, 


Hk 


= 7~Lk b 


= -1, 




M Da 


= M D = 


= M Db = 


o, 


M K 


= M Kb 


= -1/2, 





a/A 




af am Da \ 
2a a(m B — 8) 
a(2 + m,D B ) 
arriK/A 1 / 2 



ah Da /A 1 / 2 
aho/A 1 / 2 
ahnJA 1 / 2 
a(2 + h K )/A 
a(l/2 + \ikb)IA olttikb 

I A 1 ' 2 j 
(4-51) 

We now need to determine the coefficients hi and 
in order to obtain a well behaved system of evolution 
equations. 

4- Integrability and hyperbolicity 

Since the D's arise as spatial derivatives of the lapse 
and metric components, their evolution equations are ob- 
tained by taking a time derivative of their definition and 
then changing the order of the partial derivatives. If one 
later adds multiples of the hamiltonian and momentum 
constraints to the evolution equations for the D's, one 
finds that whenever these constraints are violated the D's 
in fact cease to be derivatives of metric functions. One 
consequence of this is the fact that in the sources of the 
evolution equations for the eigenfields Wi, the coefficients 
cm will no longer be invariant under rescalings of the 
form Wi — £l(a, A, B)wi (the proof presented at the end 
of section Sec. Ill C"2l of the invariance of these coefficients 
under such rescalings relied on the derivative constraints 
being satisfied). Such a property of our system of equa- 
tions is undesirable, as it makes the source criteria for 
avoiding blow-ups impossible to apply in practice. 

This leads us to the "integrability criteria", which 
states that the D's should remain derivatives of the met- 
ric functions independently of the constraints, and im- 
plies that we must set the /id's and m B s to zero and 
consider only adjustments in the evolution equations of 
the extrinsic curvature variables. 

Doing this we obtain for the characteristic ma- 
trix H4.51fl the following eigenvalues 



Ac = 0, 



±a\ 



2VA 



(m KB ± ^4: + 8h KB +Wk b ) 



(4.52) 
(4.53) 
(4.54) 
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The first three eigenvalues (Ao, A£) are precisely the ones 
we found for the 1+1 dimensional case, so clearly a£ are 
again gauge speeds (they depend on the gauge function 
/). The last two eigenvalues depend on the choices we 
have made to add constraints to the evolution equations, 
so we will call them "constraint speeds" . In fact, it is not 
surprising that we find only characteristic speeds associ- 
ated with the gauge and the constraints, since in spherical 
symmetry it is well known that there are no gravitational 
waves, i.e. no physical waves propagating at the speed of 
light. 

From the last expressions we see that if we want the 
constraint speeds to be centered along the time lines we 
must ask for 

m KB = . (4.55) 
If we now rewrite the parameter \ik b as 



2 2 7 ' 



(4.56) 



with /x a constant, then the eigenvalues A± take the fol- 
lowing simple form 



(4.57) 



At this point, only the adjustments to the evolution 
equation of the trace of the extrinsic curvature, hx and 
rriK, and the constant /x remain as free parameters. It 
is now not difficult to show that the system of evolu- 
tion equations will be strongly hyperbolic, i.e. the ma- 
trix H4.51fl can be diagonalized and a complete set of 
eigenvectors exists, as long as / > and \i £ {0, ±1}. We 
want to point out here that the latter condition implies 
that all characteristic speeds have to be distinct. How- 
ever, there is one important exception, since for \fx\ — 1 
the adjustments hx = —2 together with tjik = also 
yield a strongly hyperbolic system (which will be of some 
importance later). Furthermore we observe that no par- 
ticular problems arise for the case of harmonic slicing 
(/ = 1). For completeness we explicitly state here the 
general form of the eigenfields: 



w = D a - f -D-2fD B 



(4.58) 



w± 



v?(i 



D„± 



± 



■w. 



2(2 + h K )±m K \ff 
Mv 7 / D B ± 2VJ Kb . 



fAK 

2 



D B 

AK b , (4.59) 
(4.60) 



In the case when = 1, hx = —2 and ttik = 0, the 
eigenfields tUj. and vu± take a different form given by 



-£ = y/fAK±D a , 

= v7 d b ±2Vak b 



(4.61) 
(4.62) 



5. Indirect Linear Degeneracy 

Just as before, the linear degeneracy criteria is triv- 
ially satisfied. Applying the indirect linear degeneracy 
criteria, dXi/dwi — 0, to the eigenvalues A 7 



\£, we find 



dw± 



2(1- n*)VJA 



J 2 



(4.63) 



This gives us the same condition on / for avoiding a blow- 
up as before, namely condition (|4.17(l . which is precisely 
the result we obtained for the 1+1 dimensional case. In 
addition, however, we note that blow-ups can also arise 
from the second pair of eigenvalues, for which we find 



Owl 



J 2 



4 + 2h K ± nm K y/J 



2/ 



1- M 2 

The condition for avoiding these blow-ups is then 



(4.64) 



J 2 



4 + 2h K ±fim K y/f 



2/ 



= . 



1 -n 2 

(4.65) 

The first thing to notice is that we clearly must take 



m K = , 

in which case the condition reduces to 
'4 + 2h K 



/' 



2 



l-/x 2 



2/ 







(4.66) 



(4.67) 



Now, if we insert a member of the gauge shock avoiding 
family into this condition, we find 



M/ = 0, 



(4.68) 



which, remembering that for strong hyperbolicity one 
must have / > and /i / 0, brings us to the rather 
discouraging result that - for the adjustments considered 
here - we can not avoid both gauge shocks and constraint 
shocks coming from the indirect linear degeneracy crite- 
ria at the same time. 

Instead of using a full blown member of the gauge 
shock avoiding family, we could be less ambitious and use 
a solution that avoids gauge shocks only approximately. 
For example, in Ref. |27| it was shown that the stan- 
dard 1+log slicing, corresponding to the choice / = 2/ a, 
avoids gauge shocks to first order (which explains why it 
is so robust in practice) . Taking such a form of / we find 
that the condition for avoiding constraint shocks is 



h K = -2 



a — fx 



(4.69) 



For \/j,\ « 1 this condition implies hj< ~ —2 which, as 
we have seen in the previous section, is the only value 
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of Hk that yields a strongly hyperbolic system when 
\fi\ = 1. However, as mentioned before, for the special 
case {|/z| = 1, Iik — — 2} the eigenfields are different and 
the analysis should be performed separately. When we 
apply the indirect linear degeneracy criteria to the eigen- 
fields w± given by 1)4.62(1 . we find that 9A£/9w£ is non- 
zero, and since in this case we have no free parameters, it 
is hard to say if the condition is "closely satisfied" or not. 
Still, Eq. (|4.69[) does seem to indicate that the combina- 
tion {|//| = l,hpc = —2} is preferred by indirect linear 
degeneracy. We will return to this case when we consider 
some numerical examples below. 



6. Source Criteria 

The source criteria for the gauge eigenfields wi^ yields 
the same condition on / as indirect linear degeneracy 
since the quadratic coefficient turns out to be 



-±±± — 



±- 



2(1 



J 2 



(4.70) 



On the other hand, for the constraint eigenfields w± we 
find the following quadratic coefficient 



c ±±± 



± 



(l + /i 2 ./> 
l6v 2 fVA 

7 + 4h K -3(j, 2 f±2fj,m K y/f] . (4.71) 



If we want to avoid a blow-up this coefficient must vanish, 

7 + Ah K -3(i 2 f ±2fim K y/f = . (4.72) 
This can be accomplished for any / if we choose 



m K = , 



and 



l + 3^ 2 / 



(4.73) 



(4.74) 



We first notice that again we obtain the condition 
rriK = 0, just as we found in the previous section. It 
is also interesting to notice that for a given choice of / 
we have a one parameter family of solutions that avoids 
these type of shocks. From (|4.56l) and (|4.74|) we see that 
the parameter fi relates hx and hx B according to 



h K = -1 



3/i 



K B 



(4.75) 



Considering the restrictions on fj, imposed by strong hy- 
perbolicity, fi ^ {0,±1}, we see that in the (h,K,hK B ) 
plane this shock avoiding family must be such that 



hx B > ~ , h Ks ^ -(/ - 1) 



(4.76) 



We see hence that by appropriate choices of f(a) and 
suitable adjustments to the evolution equations, it is pos- 
sible to avoid at the same time the gauge and constraint 
shocks identified by the source criteria. 

As a final comment it is important to point out that, in 
contrast to what we found for the gauge shocks, in this 
case indirect linear degeneracy and the source criteria 
yield different statements for avoiding constraint shocks. 



C. Numerical tests in spherical symmetry 

We will now describe some numerical experiments de- 
signed to test the shock avoiding conditions found in the 
previous sections. We will concentrate on two different 
types of tests: The robust stability test j^jji and a test of 
Minkowski initial data with a violation of the constraints 
added by hand. 



1. Robust stability test 

As a first numerical experiment we have performed the 
robust stability test as described in [4p|. For this test 
one takes Minkowski initial data and adds random noise 
with a small amplitude to all dynamical variables. For 
the evolution we have used both harmonic slicing with 
/ = 1 and standard 1+log slicing with / = 2/ a. 

For the simulations discussed below we used 1,000 grid 
points, a grid spacing of Ar = 0.001 and a time step 
At = Ar/2. Furthermore, we have demanded periodic 
boundary conditions and have set all 1/r and 1 /r 2 terms 
to zero by hand (i.e. we are performing the run "at in- 
finity"). Setting these terms to zero should have no im- 
portant effect on the presence of blow-ups since one can 
readily verify that such terms are only linear in the eigen- 
fields. Moreover, one could easily regularize the equa- 
tions at the origin by following the procedure described 
in 39], but doing so would only obscure the study of 
blow-ups by mixing them with purely geometric effects. 

We start the simulations with Minkowski initial data 
plus random noise of order 10 -6 on all evolution vari- 
ables. At latter times we compute the error 5 as the 
average absolute value over the grid of the quantity 

(ELk-il + ELikD/s- 

We first performed runs for the ADM system without 
adding constraints, and observed the well known growth 
in the average error caused by the fact that ADM is only 
weakly hyperbolic. Next, we implemented the /i-family 
given by (|4.56|l and l|4.74|) , which according to the source 
criteria is shock-avoiding. Figure [5] shows the behavior 
of the average error for this family. In the top panel 
we plot the growth of the error for the case of harmonic 
slicing as a function of time (measured in units of the 
light-crossing time of our computational grid), for sev- 
eral different values of We see that for ADM (which 
is not a member of the /i- family), and for the cases with 
l-i = and \fi\ = 1, the error grows linearly with time, 
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FIG. 5: Top. Using harmonic slicing, the average error for the 
robust stability test is shown for ADM and different members 
of the /i-family as a function of time (measured in units of 
the light crossing time of our grid). Bottom. The value of 
the average error after one light crossing time is plotted as a 
function of the parameter fi. 
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FIG. 6: Top: Setting \n\ — 1, the average error for har- 
monic slicing is found to grow only moderately for Ilk = —2, 
whereas for other values of Hk (and in particular for ADM 
corresponding to Kk — 0) a rapid growth is observed. Bot- 
tom: The average error after one crossing time is plotted as 
a function of ha- 



while for \n\ = 1/2 and = 2 the error initially does not 
grow at all (at later times, however, also in these cases a 
linear growth with a very small gradient develops). The 
lower panel in this figure shows the value of the average 
error after one light-crossing time as a function of |//|, as 
obtained both for harmonic and l+log slicings. We see 
that for values of close to or 1, the error is already 
very large after one light crossing time, while away from 
these values the error remains small. The poor behav- 
ior of the simulations with fi — and = 1 is caused 
by the fact that for such cases the evolution system is 
not strongly hyperbolic. Also, from (|4.57|) we see that 
for values of close to but different from either or 1, 
the eigenspeeds \± associated with the constraint modes 
become very similar to either Ao or Aj_ . This means that 
even if the system is still strongly hyperbolic, the argu- 
ment used for ignoring mixed terms in the sources will 
not apply. 

In particular, for the case \fi\ = 1 the adjustment 
1%k = — 2 suggested by the indirect linear degeneracy cri- 
teria turned out to be helpful. Because of this in Figure^] 



we also show similar plots testing different values of the 
parameter hx in the case of |ju| = 1. Here one can observe 
that for values of \ik other than -2, and in particular for 
ADM corresponding to = 0, a linear growth in the 
average error is present. For the adjustment suggested by 
the indirect linear degeneracy criteria (which is the only 
value here which yields a strongly hyperbolic system), 
initially no error growth is found. 



2. Minkowski initial data plus constraint violation 

As we saw in the previous section, the robust stabil- 
ity test is very good at distinguishing between strongly 
and weakly hyperbolic systems, but does not show any 
clear indication that, among strongly hyperbolic systems, 
those that avoid shocks are better behaved. This should 
not be surprising as the robust stability test uses ran- 
dom, uncorrelated and non-smooth initial data, while 
shock formation is the result of the coherent evolution 
of smooth initial data. 
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For this reason, using harmonic slicing we have per- 
formed evolutions of Minkowski initial data with a 
smooth violation of the constraints added by hand. Here 
we have chosen a perturbation in Kb similar to the one 
we used for the scalar wave equation, namely 

K B {t = 0) = - (2nr/a 2 ) exp {-r 2 / a 2 ) , (4.77) 

with parameters k — 0.05 and a = 1. 

For the simulations shown below we used 8,000 grid 
points and a grid spacing of Ar = 0.01 (which places the 
boundaries at ±40) together with a time step At = Ar/2. 
Furthermore, we have again neglected 1 jr and 1/r 2 terms 
as the simulation can be "shifted" to large values of r. 

In Figure [7] we show contour plots of the root mean 
square (r.m.s.) of the Hamiltonian constraint as a func- 
tion of the adjustment parameters and hx B at two 
different times during the evolution, using 40 equidistant 
parameter choices in each direction. The momentum con- 
straint, not shown here, has a very similar behavior. Note 
that black and dark gray colors correspond to regions 
where the r.m.s. of the Hamiltonian constraint grows 
rapidly, and light gray denotes regions where it grows 
very slowly or not at all. For ADM (corresponding to 
hjc = hx B = and denoted by a black circle) a rapid 
growth of the constraints is found. We also observe rapid 
growth close to the white circle, representing the special 
case with hx — —2 and h,K B = which corresponds to 
the only strongly hyperbolic system along the line = 1 
and is preferred by indirect linear degeneracy. 

The white diagonal line corresponds to the shock 
avoiding ^-family obtained from the source criteria (not 
defined for the points /i = and = 1 represented as 
boxes). We clearly see that this one-parameter family 
falls in the middle of the region where the r.m.s. of the 
Hamiltonian constraint does not grow, indicating that it 
does correspond to a preferred region of parameter space. 

There is a final point related to Figure The figure 
shows that the line = 1 also seems to produce slow 
growth of the constraints. However, as mentioned before, 
in that case the system is only weakly hyperbolic and our 
whole analysis breaks down, so we have no explanation 
as to why this line represents a preferred region. 

In order to see the formation of shocks more clearly, 
in Figure |S] we show the time evolution (shown every 
At = 2) of the eigenfields w± which are associated with 
the formation of constraint shocks. The upper panel 
shows the evolution for the parameter choice hx = — 2 
and h,K B = 1/2 (which for harmonic slicing implies 
= \/2), corresponding to a strongly hyperbolic case 
that does not avoid constraint shocks. From the figure we 
see how a shock is clearly forming, as expected. The mid- 
dle panel corresponds to the parameters Hk — —2 and 
hx B = (and hence = 1) preferred by the indirect lin- 
ear degeneracy criteria. Again we observe the formation 
of shocks. Finally, the lower panel shows the evolution 
for a member of the shock-avoiding /i-family with param- 
eters hx = —1/4 and hx B —1/2 (i.e., |^| = V2). In this 
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FIG. 7: Contour plots of the r.m.s. of the Hamiltonian con- 
straint are shown at times t = 6 and t = 18. 



case the evolution shows a wave-like character, and no 
shocks form on the time-scale considered here. 



V. CONCLUSIONS 

We have presented an analysis of two different blow- 
up mechanisms for systems of hyperbolic equations of 
the type found in general relativity, namely the geo- 
metric blow-up or gradient catastrophe, and the ODE- 
mechanism. As an example of how these mechanisms 
operate we have used the simple one-dimensional wave 
equation with dynamic wave speed and non-trivial source 
terms. 

We have later performed a blow-up analysis of one- 
dimensional systems in general relativity, concentrating 
on "toy" 1+1 relativity and spherically symmetric rel- 
ativity, and using the hyperbolic Bona-Masso family of 
slicing conditions. In the first case we have recovered the 
well known gauge shocks and found, somewhat surpris- 
ingly, that both blow-up criteria give precisely the same 
condition for shock avoidance. When studying the spher- 
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FIG. 8: Evolution of the eigenfields (black) and w'L (gray) 
associated with constraint shocks. Top. Blow-ups occur for a 
strongly hyperbolic case that is not shock avoiding. Middle. 
Shocks also form for the adjustments suggested by the indirect 
linear degeneracy criteria. Bottom. No shocks are found for 
a member of the /i-family obtained by the source criteria. 



nally, we have presented numerical simulations that con- 
firm that constraint shocks can and do indeed form, and 
that they can be avoided by using the results of our pre- 
vious analysis. 

We would like to mention here that although in the 
present study we did not work with adjustments that 
can avoid constraint shocks according to both blow-up 
criteria, this is simply because we have considered only a 
restricted form of adjustments to the evolution equations. 
More elaborate adjustments can give rise to completely 
shock free formulations [45| . Moreover, empirically we 
have found that the source criteria seems to be far more 
important than indirect linear degeneracy (though it is 
still not clear why in the case of gauge shocks both crite- 
ria give rise to the same condition for avoiding blow-ups) . 

As a final note we want to point out that there is a 
very important difference between gauge and constraint 
shocks. When a gauge shock develops, it means that 
our coordinate system has broken down and there is 
no way to continue the evolution past that point other 
than choosing a different gauge. Gauge shocks are cer- 
tainly non-physical, as the geometry of spacetime re- 
mains smooth, but they do represent very real patholo- 
gies in the spatial foliation. Constraint shocks, on the 
other hand, are not only non-physical, but they should 
not arise at all if the constraints remain satisfied exactly. 
In a numerical simulation the constraints are of course 
violated, but as the resolution is increased this violation 
should become smaller and smaller, so the possible forma- 
tion of constraint shocks becomes less of an issue. Nev- 
ertheless, at any fixed resolution constraint shocks can 
form at a finite time unless care is taken to use a form of 
the equations that avoids them. 

Because of this we recommend that out of the many 
possible ways of constructing strongly hyperbolic evolu- 
tion systems in general relativity, for numerical evolu- 
tions one should concentrate on those that have better 
shock avoiding properties, both in the gauge and con- 
straint violating sectors, as these should prove to be more 
robust in practice. 
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